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Abstract. In this paper we study the realization of lattice models in mixtures 
of atomic and dipolar molecular quantum gases. We consider a situation where 
polar molecules form a self-assembled dipolar lattice, in which atoms or molecules 
of a second species can move and scatter. We describe the system dynamics in 
a master equation approach in the Brownian motion limit of slow particles and 
fast phonons, which we find appropriate for our system. In a wide regime of 
parameters, the reduced dynamics of the particles leads to physical realizations 
of extended Hubbard models with tuneable long-range interactions mediated by 
crystal phonons. This extends the notion of quantum simulation of strongly 
correlated systems with cold atoms and molecules to include phonon-dynamics, 
where all coupling parameters can be controlled by external fields. 
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1. Introduction 



The recent preparation of cold ensembles of homonuclear and heteronuclear molecules 
in the electronic and vibrational ground state 




EQEBESBBSSEiH S3 Bill; El lUlii [IlTii has opened 

the door to a new chapter in the theoretical and exp erimental study of trap ped 



guantum degenerate g^esJ3llj32j; l33|; [3^; l35|; l36t l37|; l38|; l39|; l40|; l41|; [4^; [43|; [44 



3 4(| 47; 48; 4^; 5(| 51; 53" 53| . Heteronuclear molecules, in particular, have large 
electric dipole moments associated with rotational excitations, and the new aspect 
in quantum gases of cold polar molecules is thus the large, anisotropic dipole-dipole 
interactions between the molecules which can be manipulated via external DC and 
AC fields in the microwave regime 31; 32; 33|; 34; 52|; 53; 54 1. In combination with 
reduced trapping geometries, this promises the realization of novel quantum phases 
and quantum phase transitions, for example in the case of bosons the transition from 
a dipolar superfluid to a crystalline regime [13] ■ The theory of dipolar quantum gases, 
and various aspects of strongly correlated systems of polar molecules has been recently 
reviewed by Baranov [55| and by Pupillo et al. [5|| 

Below we will extend these theoretical studies to mixtures of atomic and dipolar 
molecular quantum gases. More specifically, we will be interested in a situation 
where polar molecules form a self-assembled dipolar lattice, in which atoms can move 
and scatter. This scenario leads to a new physical realization of a Hubbard model 
where atoms see the periodic structure provided by the crystal formed by the polar 
molecules. In contrast to the familiar case of the realization of Hubbard models 
with cold atoms in optical lattices, where standing laser light waves produce a fixed 
periodic external lattice potential, self-assembled dipolar lattices have their own lattice 
dynamics represented by phonons. Thus atoms moving in dipolar crystals give rise to 
Hubbard models which include both (i) atom-phonon couplings, and (ii) atom-atom 
interactions. This extends the notion of quantum simulation of strongly correlated 
systems with cold atoms to include phonon-dynamics, where both the atomic and 
phonon coupling parameters can be controlled by external fields. We note that atom - 
molecule mixtures arise naturally in photoassociation experiments, where a two-species 
atomic quantum degenerate gas is partially transferred to ground state molecules via 
formation of highlying Feshbach molecules, followed by a Raman transfer to the ground 
state. Similar Hubbard models result also in mixtures of two unbalanced species of 
polar molecules, where the first molecular species forms a crystal while the second 
species provides the extra particles hopping in the self-assembled dipolar lattice. 

In a recent work [571 ] , we have shown that the dynamics of atoms and molecules 
embedded in dipolar crystals is conveniently described by a polaronic picture where 
the particles are dressed by the crystal phonons. Standard treatments of polaron dy- 
namics show a competition between coherent and incoherent hopping of a particle in 
the lattice. The former corresponds to tunneling of a particle from one site to to the 
next, while carrying the lattice distortion, without changing the phonon occupation. 
The latter corresponds to thermally activated particle hopping, related to incoherent 
hopping events where the number of phonons changes in the hop. That is, the polaron 
loses its phase coherence via the emission or absorption of phonons [5a ] . The physics 
of polarons has a long history, dating back to the seminal work of Landau I59H . Ex- 




cellent reviews on the subject can be found e.g. in [6C 
Here we take the simple approach of placing these coherent and incoherent processes 
in the natural framework of a master equation treatment of the system dynamics, 
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Figure 1. A dipolar crystal of polar molecules in 2D (a) and ID (b,c) provides 
a periodic lattice V cp for extra atoms or molecules giving rise to a lattice model 
with hopping J and long-range interactions Vi j (see text), (a) In 2D a triangular 
lattice is formed by polar molecules with dipole moment d c perpendicular to the 
plane. A second molecular species with dipole moment d p -C d c moves in the 
honeycomb lattice V cp (darker shading corresponds to deeper potentials), (b) A 
ID setup with atoms scattering form a dipolar crystal with lattice spacing a. (c) 
A ID dipolar crystal provides a periodic potential for a second molecular species 
moving in a parallel tube at distance b. 



where the phonons are treated as a thermal heat bath. We extend our work [57| by 
re-deriving the results of [58| in the master equation context, and calculating addi- 
tional corrections to the coherent and incoherent time evolution for our atomic and 
molecular mixtures. Since we find that in the latter the polaron dynamics is typically 
slow compared to the characteristic time evolution of the bath, we specialize to the 
Brownian motion limit of the master equation |68|; |69| . We show that for the models 
of interest and low-enough temperatures, corrections to the coherent time evolution 
of the polaron system are small, and thus the dynamics of the dressed particles is 
well described by an effective extended Hubbard model in a wide range of realistic 
parameters. 

The paper is organized as follows. In Section [2] below we derive the Brownian 
motion master equation for our system, providing explicit expressions for the 
coherent and incoherent polaron dynamics (details of the derivation can be found 
in |Appcndix C[ ) . In Section [3] we review the realization of dipolar crystals in two and 
one dimensions. In Section [4] we provide details of the implementation of polaronic 
extended Hubbard models for two one-dimensional configurations with atoms and 
molecules. 

2. Dynamics of particles trapped in a crystal of polar molecules 

We show below in Section[4]and |Appcndix "A| that in a wide range of system parameters 
the dynamics of atoms or molecules which move in a self-assembled lattice of dipoles 
is well described by the following Hamiltonian 
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+ huj ^ a Lx a ^ + M q ^e l ^c] Cj {a^ x + al q A ), (1) 




where the brackets () indicate sums over nearest neighbors. The latter is a single- 
band Hubbard model for the particles coupled to the acoustic phonons of the lattice. 
The first two terms on the r.h.s. of ((T|) describe the hopping of particles between 
neighboring sites of the lattice with a tunneling rate J, and the density-density 
interactions with strength Vij for particles at site i and j, respectively. Here, c\ (c,) 
denotes the creation (annihilation) operator for a particle at site i. These particles 
can be either fermions or bosons. The third term describes the excitations of the 
crystal given by acoustic phonons, where A (aq.A) creates (destroys) a phonon with 
quasimomentum q in the mode A, with dispersion relation LUq_\. The last term in 
([I]) is the particle-phonon coupling, which is of the density-displacement type, with 
coupling constant Mq.\. The microscopic derivation of ([I]) is detailed in | Appendix A| 
and | Appendix B| For the models we consider (see Section [3] and Section |4] below), we 



(i) The phonon coupling can largely exceed the hopping rate J, which precludes a 
naive treatment of the particle-phonon coupling as a (small) perturbation; 

(ii) We are generally interested in the so-called non-adiabatic regime, where the 
characteristic phonon frequency Hlod (the Debye frequency) is typically (much) 
larger than the average kinetic energy of the particles ~ J, that is Hlod 3> J, (see 
Section [4]). This is due to the fact that, unlike e.g. the case of electrons in ionic 
crystals, in our system the mass of the particles and of the crystal dipoles are 
comparable, and the crystals can be made stiff. 

Below, we derive a master equation for the dynamics of the particles only, while 
the crystal phonons are treated as a thermal heat bath. The master equation in the 
Markovian limit has the general form 



where H$ denotes a Hamiltonian term for the coherent time evolution of the reduced 
density matrix p§ for the particles, while T> describes dissipative processes responsible 
for the incoherent dynamics. The coherent processes correspond to hopping of a 
particle dressed by the crystal phonons (a polaron) and polaron-polaron interactions, 
while the latter are thermally activated incoherent hopping, where a particle loses 
its phase coherence by emitting or absorbing a phonon in the hopping process. 
As suggested by points i) and ii) above, we derive J2} in a perturbative, strong- 
coupling, approach where the hopping rate J in JTJ) acts as the small parameter. In 
addition, we work in the Brownian motion limit of the master equation, where the 
characteristic time evolution of the system rg ~ max(l/J, 1/Vij) is (much) slower than 
the characteristic time evolution of the thermal heat bath tb ~ l/fuo-Q. We obtain 
explicit expressions for the coherent and incoherent (dissipative) contributions to the 
system time evolution, and show that the latter can be made negligible in a wide range 
of realistic parameters for our models. 

We find that the effective Hamiltonian for the coherent dynamics of the particles 
only has the form 



find that: 



p s (t) = -~[H s ,ps(t)] +V[ Ps {t)], 



(2) 




(3) 
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which corresponds to an extended Hubbard model. Here J and Vy are the tunneling 
rate for the particles and their mutual interactions, respectively, which are modified 
with respect to their bare values in ([T]) by the coupling to the crystal phonons. 

2.1. Lang-Firsov Transformation (Polaronic picture) 

In our system (see Section [4]) , we find that particles are slow and strongly coupled 



to the crystal phonons, with Hamiltonian ((T|). Following Lang and Firsov 6l|, it is 
convenient to change from this picture of bare particles strongly coupled to phonons 
to an equivalent one, where particles freely hop in the lattice while carrying the lattice 
distortion. This corresponds to dressin g th e particles with the lattice phonons, and 
the dressed particles are named polarons [581 ] . This is achieved by performing a unitary 
transformation of H as H = UHU^ with 

U = exp[£ .i.^'^ c rj.r^ - a q , A )] (4) 

and Uq t x = M a .\/huj a .\ displacement amplitudes. In the new picture the Hamiltonian 
reads [HlH 



+ \^ViA c \ c i c i> ( 5 ) 

where now Cj (cj) are the annihilation (creation) operators of a polaron located at site 
j. Here, 

X, = exp[- J2 Uq,Ae iqr V_ qiA ~ «q,A)]- (6) 

is the displacement operator for the crystal molecules due to the presence of a particle 
located at site j. The energy E p in |5]) is the polaron shift 58] 

£^, = X;<a/Ka (7) 

q,A 

and N p = Y) j CjCj is the total number of particles. The quantity 

M 2 

V,, = V± 4 - 2 J2 cos[q(r° - r°)] = V ij + (8) 

q nuj^x 

is a modified particle-particle interaction, which comprises two terms: The first is the 
original (bare) interaction, while the latter is the phonon mediated particle-particle 
interaction V^j, which in general (i) is long-ranged and (ii) can be comparable in 
strength to the bare interactions V%j , see Section [4] 

For J = the new Hamiltonian ([5|) is diagonal and describes interacting polarons 
and independent phonons. The latter are vibrations of the lattice molecules around 
new equilibrium positions with unchanged frequencies. For small finite J, we can treat 
the first term in (O perturbatively, which is the starting point for the master equation 
approach detailed in the following sections. 
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2. 2. Effective dynamics of polarons inside a thermal crystal 

In this section we derive a master equation describing the effective dynamics of the 
polarons embedded in the crystal, which we assume to be in thermal equilibrium. 
That is, we consider the phonons to provide a heat bath at temperature T with a 
density matrix given by 

P°B=Un^T)exp(-^al x a^Y (9) 

q.A 

Here n^\(T) denotes the (Bose-Einstein) phonon distribution at temperature T, 



ftq,A(T) = (a q , A a q ,A> = e _^ q . A/ 1 fcBT _ 1 - 



where (O) = tr^iOp^} is the thermal expectation value of the bath operator O, with 
tr B denoting the trace over the bath degrees of freedom. 

We split the Hamiltonian ([5]) into three parts as H = H$ + Hb + H\ with 

H S = - J J2 c\c 3 (Xjxj) + ViAcfa*' ( 10fl ) 
Hb = E ^q,A a q.A a q,A' ( 10& ) 

q,A 

H i = - J E c h(x}Xi - (xjx,)). (ioc) 

(i,3) 

Here, H$ is the "reduced system" Hamiltonian describing the dynamics of polarons, 
with a hopping rate J(xjXj) and interactions Vij for two polarons at site i and 
j. The Hamiltonian H# is the Hamiltonian for the bath, and Hamiltonian Hi gives 
the interactions between the system and the bath. In writing (|10aj) - (I10c|) we have 
conveniently added a term —JJ2{ij) c \ c j(X\Xj) to H$ and subtracted it from H\. 
This ensures that the thermal average over the interaction Hamiltonian vanishes, 

(i?i) = o|s3. 

The expectation value (XjXj) in equation (|10a|) reads (see |Appendix CTj ) 

q(r° -r )" 

(2n q , A (T) + 1) . (11) 



(XjXj) = e- s ' T with S T = 2 ^ w q . A sin 2 

q,A 



Equation <\10a\i shows that 

J=J(x}x j ) = Je~ s - 

plays the role of a phonon-modified tunneling rate, which is exponentially suppressed 
for St > 0. In the following we will often distinguish between two regimes, i.e. a 
weak coupling regime where St -C 1 (and J ~ J) and a strong coupling regime where 
S T > 1 (and J < J). 



2.2.1. The Quantum Brownian Motion Master Equation (QBMME) In this section 
we derive a master equation for the coherent and incoherent dynamics of the reduced 
system of interacting polarons in the Brownian motion limit, where the system time 
evolution, rg, is slow compared to the characteristic evolution time of the bath, Tq. 
That this approach may provide a reasonable description of the system is suggested by 
the fact that in a wide range of parameters for atoms and molecules trapped in dipolar 
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crystals we find rg ~ max(l/J, l/Vij) 3> tb ~ 1/wd- We will thus derive the coherent 
time evolution and provide explicit analytic expressions for the corrections both 
to the coherent and incoherent dynamics. In Section0]we show that for the models of 
interest these corrections are in fact negligible in a wide regime of parameters, which 
provides an a posteriori self-consistency check for the approximations made here. 

The time evolution for the density matrix of the entire system p(t) comprising 
the (polaronic) reduced system and the heat bath is dictated by the Liouville-von 
Neumann equation 

p(t) = -~ [Hi(t),p{t) 

where A(t) = g^ot/Hj^g-iHot/H denotes an operator A in the interaction picture with 
respect to H = Hg + Hb . 

We assume that the reduced system and the heat bath are uncoupled at time t = 0, 
so that p(t = 0) can be written as the tensor product p(0) — ps(0) ® p B , with ps(0) 
and p B the density matrices of the reduced system and the bath, respectively. The 
condition J <C fkt>x), which states that the interaction is a weak perturbation, forms 
the basis for a Born-Markov approximation with the phonons a finite temperature 
heat bath, providing the following master equation for the reduced density operator 
of the polarons 

1 C°° - 

dTtr B {[Hiit) i [H I (t-T),fa{t)®&]]} (12) 



o 



J2 poo 

v \ dT J2 (4-(^ T ) 5 l(*) 5 i(*)[4(*-^)^(*--r),psW] 

J ° (i,i),<M) 

$ (-r, T) [cl (t - r)cj (t - r), p s (*)] cj (f )c, (t)) , (13) 



where ffl (r, T) are the (complex) bath correlation functions, 

$j(r,T) = (x+(^(i)xl(tO*K*0) - (x/(t)Xi(0><^(O^(*')> |f=t-r 

= {XjX 3 Xt(~r)M~r)) - e" 2 ^ = ^(-r,T)*, (14) 

which are treated in detail in Section 12.2.21 and | Appendix C.l| 

Under the additional condition max(J, Vij) <C (/zwdj^b? 1 ), we focus on the 
Brownian motion limit of (|13|) . where the system, like the interaction, evolves on 
a timescale much slower than the bath, see e.g. & In this limit we can approximate 
the system operators in equation (|13p by 

7 

Cj(-T) « Cj - -[H s , Cj-Jr, (15) 
and the master equation takes the form 
p s(t)= _i[ J ff s _ £ Ag(T)clc,-4 Q ,psW 

- E r fj( T )(4ci4 c i/°s(*) + PsWcJcjcJq - 2c{c lPs (t)4c 3 

(i,j){k,l) 



7#(T) [c\c h [J( £ - 2 c ^') + - K lfe )4c m cLQ, pa(t) 



(i,j)(W> fe 
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- * E^( T )[ c i c ^[ J (E4'^-E4^)+E(^-^)4 c ™4Q,Psw]],(i6) 

(i,j){k,l) k' V m 

with the quantities A|- (T) , r*j (T) , <5,fJ (T) and 7§(T) given by 

t2 /-oo 

Ag(T) = — jf drlm[^(r,T)], (17a) 

t2 /-oo 

r ^ T ) =-jpj drRe[eg(T,r)], (176) 

t2 /-oo 

7$C0 = tf/ o drrRe[^(r,T)], (17c) 

t2 /-oo 

Here Re[/(x)] and Im[/(x)] denote real and imaginary part of /(x). 

The first term inside the commutator on the r.h.s. of (|16p describes the coherent 
time evolution for the reduced system, with H§ the effective Hubbard Hamiltonian 
(I10a| . The terms proportional to A^J(T) in (fT6|) are self-energies, which in the 
single-polaron limit provide both a shift to the ground-state energy, and next-nearest- 
neighbor hopping terms [7lj |. In addition, in the many-polaron problem they can 
provide offsite polaron-polaron interactions, whose strength will be evaluated in 
Section [2 . 2 . 31 below. The terms on the r.h.s. of (p~6|) proportional to Tfj(T) are related 
to incoherent hopping events where the number of phonons changes in the hop. That 
is, the polaron loses its phase coherence via the emission or absorption of phonons. 
These processes are thermally activated and can dominate over the coherent hopping 
rate J for large enough temperatures T (58|. 

The terms proportional to TyfT) and S£j(T) are (small) corrections to the 
coherent and incoherent time evolution, respectively, which derive from the term 



proportional to r in the expansion (|15|) . and thus correspond to the Brownian motion 
corrections to the time evolution of the reduced system. 

We notice that the terms proportional to (|17a[) - (|17ril) in (JTHJ) do not identify 
directly the corrections to the coherent time evolution given by (|10a[) , because 
equation (|16p is not diagonal. Instead these terms are elements of matrices whose 
eigenvalues are the corrections. The diagonalization of the master equation can be 
done analytically in the single polaron limit and is shown below in Section 12.2.41 

In Section 12.2.31 below we provide analytic expressions for the coherent and 
incoherent corrections to the coherent time evolution given by H$, cf. (|T0a|) . In 
Section [4] we show that these corrections are in fact negligible in a wide range of 
realistic parameters for atoms and molecules, and thus H$ properly describes the 
dynamics of polarons inside the dipolar crystal. 



2.2.2. Correlation Functions: The bath correlation functions (t, 7 1 ) appearing in 
the master equation (|16[) are computed in | Appendix C.i] and read 



£«( r ,T) = e- 2S -(e-*^' T >-l), 



with *?At,T) 



dwjM(w) 



coth 



/ hw 



\k B T 



cos(wt) — isin(wr) 



(18) 
(19) 
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Here we have introduced the spectral density 

j%(w) = v bz J2 [ rfq^ 1 

A J 

where Vbz is the volume of the Brillouin zone, q d ~ 1 denotes all components of the 
quasi- momentum vector except q x , and gfj (q) reads 

ffy(q)= co S [qK -r fe |]-cos[q|rO-r°|] 

-cos[q|r°-r°|]+cos[q|r°-r°|]. (21) 

The quantity is a decaying function of the time r with max(|$^-|) = 2St, the 
actual decay rate depending strongly on the spectral density of the model [58j. In 
Sections 12.2.31 below we show that at small finite temperatures this decay rate is 
fast enough to ensure that the corrections to the coherent time evolution in the 
QBMME (|16[) for our ID models are finite and small. This provides for an a posteriori 
check of the applicability of the QBMME to the polaron problem. 



du, 



q,A 



(20) 



2.2.3. Self energies and dissipation rates in the strong and weak coupling limits In 
this subsection we provide analytical approximate expressions for the matrix elements 
Ay, r*-, 7$ and 8™ in the limits of strong and weak coupling St >■ 1 and SfCl, 
respectively. The details of the performed approximations are given in |Appendix C 
while explicit results for our ID models are shown below in Section 2] Here we 
concentrate in particular on the leading self-energy corrections to the coherent time 
evolution determined by H$ in (jlOal) in the two regimes. 

In accordance with literature [6l| . we find that in the strong coupling limit, 
> 1, the self-energies are suppressed by a factor of the order of (J/E p ) 2 , with 
E p the polaron shift ([?]). In addition we find that in the weak coupling limit, St "C 1, 
these corrections are strongly suppressed by a factor (J/htur>) 2 and as a consequence, 
in Section U we show that they are negligible in a wide range of realistic parameters 
for our models. 



It is shown in |Appcndix C| that all matrix elements A^j , T 



7y and 5™ in the 

strong coupling limit are stron gly suppressed by an exponential factor of the order 
e~ 2S ' T , unless i = I and j = k, 58|; |6lJ. This makes sense, since processes for which 
i =/= I and j =/= k correspond to a double hop of a polaron, each one being suppressed 
by the factor e~ s ' T . The conditions i = I and j = k describe a "swap" -process, which 
corresponds, e.g. in the case of A^J, to the virtual hop of a polaron from its current 
position to a neighboring site and back. 
The various corrections read 



(22a) 
(226) 
(22c) 
(22 d) 




n l/2 Be -B 74Ar Erfi ( B /2V3?) 



TO 



2A T 4 
J 2 ttV2 Be-^l^ 



3/2 



h 2 ojl 



A 



3/2 
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where Erfi denotes the error function and where we have introduced the quantities 
At = J dw l -J™{w)w 2 coth (j^j , (23) 

B = J dwJ$(w)w. (24) 

Equations (|22a|) - (|22dp result from an expansion of the function — <5>^j(r, T) appearing 
in the exponent of £y (t, T) , cf. §TE\i , up to second order in the time r around its 
maximum. For a vanishing At, which corresponds to neglecting second order terms 
in the expansion of — 3>S(t, T), we find 

A - - (25) 
* = (S^F <26 » 



A simple estimate of B, see | Appendix C.2 shows that it is of the order of B ~ 
2E p /hujB for a sufficiently strong coupling and therefore Aq5/-E p cx (J/E p ) 2 and 
792* cx (J/Ep) 2 . This dependence is known in literature as the "1/A" strong coupling 
expansion, where A = E p /J, see |6lj |. 

In the weak coupling limit, St *C 1, the correlation functions £y can be expanded 
to first order in (r, T) , which leads to the following expressions for the matrix 
elements 

A «(T) » ^- e - 2S -p[ dw^l (27a) 
13 y hw D J w 

W Or.) 



as detailed in Appendix C.2 Here, dx denotes the Cauchy principal value integral. 
In contrast to the strong coupling regime, as shown in |Appendix C| here the (small) 
parameter that gives the approximate size of the corrections is J/htun and not J/E p . 
Explicit results for the corrections in this limit are given in Section 0] below. The ratio 
(J /huii)) 2 also determines the size of the incoherent processes HTjl . This can be seen 
by performing a low temperature approximation of expression (|27fc[) above, for which 
we find HT^j(T) cx J 2 k-^T / (hwo) 2 . We notice that, in addition to being proportional 
to the small factor (J /hw-p) 2 , these corrections depend linearly on temperature, a 
result also found in [58|. Because dipolar crystals can have a large Debye frequency, 
in our models we find J/hwi, <C 1 and thus corrections to the coherent time evolution 
determined by H$ are small. 



2.2.4- Diagonalization of the master equation As noted above the corrections to the 
master equation rfJ(T), A^(T),j^ l (T) and 5^(T) are just matrix elements, while the 
actual energies and rates are obtained by diagonalizing the various terms in (| 16[) . 
In the following we sketch how to perform this diagonalization in the simple case 
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of a single polaron for A*j- and 1^, while similar computations for all coherent and 
incoherent corrections are shown in | Appendix C.3| In the case of a single polaron 
the terms proportional to A*- are easily diagonalized, since the eigenergies for a 
particle in a lattice are readily determined by the energies associated with the various 
quasimomenta. That is 

£ AgCThMc =EE A^ +m+ "4c l+m+ „ (28) 

(ij)(kl) i m,n 

= ^^A™ m m+ "e^^) CqCq (29) 

q m,n 

where the sums over m,n range over basis vectors in the lattice, and y7 c]cj+ m = 
yj q e lqr ">c q Cq. The eigenvalues can now be directly read-off from PS)) , as A q = 

E m n A™^ + V^ r °«>+ r °). The largest eigenvalue gives an upper bound to the energy 
of this term in the master equation. 

To determine the rate of the dissipative term in (|16[) we notice that this term can 
be written as 

~ HtiPK^cl^AWJ-^ftC)^)' (3°) 

(ij)(kl) 

The amplitude of the dominant rate can be now estimated from Yluj){kl) ^ij CO c l c j c l c *! 
whose eigenvalues read 

r q = £ r™r + "e 4q(r '" +r " ) - (31) 

m,n 

Calculations similar to the ones above lead to the eigenvalues 7 q and <5 q for the 
Brownian motion corrections (see | Appendix C.3 ) . 

In one dimension, which is relevant for the models of Section [4] below, we find 

(32) 

(33) 

7o 1 r 1 (T))cos( 9a )], (34) 
<5 V 1 (T))cosM]. (35) 



In conclusion, we note that the corrections to the coherent time evolution given 
by A*-, rfj, J® and 5^ can be made small both in the strong and weak coupling 
regimes, by ensuring that the ratios J/E p and J/TujJt> are small, respectively. The 
smallness of these corrections provides for an a posteriori check of the applicability of 
the Brownian motion master equation approach to the polaron problem. 

3. Crystals of polar molecules 

In this section we briefl y r eview how to realize self-assembled crystals of polar 
molecules. Following [52j; l57|, here we focus on crystals in two and one dimensions. 
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(a) 




normal 




Figure 2. (a) System setup: Polar molecules are trapped in the (x, y)-plane by 
an optical lattice made of two counter-propagating laser beams with wavevectors 
±1cl = ±fcL e zi (blue arrows). The dipoles are aligned in the z-direction by a DC 
electric field Euc = ^DC e z (red arrow), (b) Phase diagram in the T — plane: 
crystalline phase for interactions > r™ and temperatures below the classical 
melting temperature T m (dashed line) 17211 , The crossover to the unstable regime 
for small replusion and finite confinement luj_ is indicated (hatched region). 



However, self-assembled crystals in three-dimensions can also be realized as detailed 
in [H3. 

Two-dimensional crystals: We consider a system of cold polar molecules in the 
presence of a DC electric field under strong transverse confinement, as illustrated in 
figure HJa). A weak DC field along the z-direction induces a dipole moment d c in 
the ground state of each molecule. Thus, the molecules interact via the dipole-dipole 
interaction (r) = D(r 2 — 3z 2 )/r 5 , with D = d 2 . This interaction is purely repulsive 
for molecules confined to the x,y-plane, while it is attractive for z > r/^/3, leading 
to an instability towards collapse in the many body system. In reference [52f it is 
shown that this instability can be suppressed for a sufficiently strong 2D confinement 
along z, as provided, for example, by a deep optical lattice with frequency fku± (blue 
arrows in the figure). In fact, a strong confinement with > D/a 3 , with 
mean interparticle distance, confines the molecules to distances larger than 

1/5 



the 



12£_ 

k m c uj±_ 

where V^(r) is purely repulsive 
m c is the mass of the molecules, 
described by the Hamiltonian 



(36) 

and thus the system is collisionally stable. Here, 
The 2D dynamics in this pancake configuration is 



tt2D 



P 2 



(37) 



which is obtained by integrating out the fast z-motion. Equation Q37p is the sum of 
the 2D kinetic energy in the a;,y-planc and an effective repulsive 2D dipolar interaction 

V™{p)=D/p\ (38) 

with pij = (xj — Xi,yj — Hi) a vector in the cc,y-plane. Tuning the induced dipole 
moment d c drives the system from a weakly interacting gas (a 2D superfluid in the 
case of bosons), to a crystalline phase in the limit of strong repulsive dipole-dipole 
interactions. This crystalline phase corresponds to the limit of strong repulsion where 
particles undergo small oscillations around their equilibrium positions. The strength 
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of the interactions is characterized by the ratio r^ of the interaction energy over the 
kinetic energy at the mean interparticle distance a 

r a - E ^ - D / a3 - Dm . (39) 
-E-kin H 2 /ma 2 h 2 a 

This parameter is tunable as a function of d c from small to large r^. A crystal forms 

for 

r d > r c = 18 ± 4, (40) 

where the interactions are dominant [52I : [73[ . For a dipolar crystal, this is the limit of 
large densities, as opposed to Wigner crystals with 1/r- Coulomb interactions. 

Figure [2b) shows a schematic phase diagram for a dipolar gas of bosonic 
molecules in 2D as a function of and temperature T. In the limit of strong in- 
teractions r^ > r c the polar molecules are in a crystalline phase for temperatures 
T < T m with T m w 0.09D/a 3 ~ 0.018r d ^ R , c , with E R>C = ir 2 h 2 /2ma 2 [z| the crystal 
recoil energy, typically a few to tens of kHz. The configuration with minimal energy is 
a triangular lattice with spacing a L = (4/3) 1/,4 a see [52j. Excitations of the crystal are 
acoustic phonons with Hamiltonian given by equation (|10b[) , and characteristic Debye 
frequency Hlod ~ l.&^fr2E^, c . The dispersion relation for the phonon excitations is 
obtained in |Appendix B.l."2| and shown in figure 0(b) below. 

One- dimensional crystals: One dimensional analogues of the 2D crystals can be 
realized by ad ding an additional in-plane optical confinement to the configuration of 
figure [^ a) [57; 74; 75|. For large enough interactions r c i ^> 1, the phonon frequencies 



have the simple form hw q = (2/tt 2 ) [I2r d f q ] 1/2 # R)C , with f q = J2j>o 4 sin (W/2) 2 /j 5 , 
see figure G2a) . The Debye frequency is huj£> = fiw„/ a ~ lAy^E^^, while the clas- 
sical melting temperature can be estimated to be of the order of T m ~ 0.2rdEi{ tC /kB, 
see [zl. 

Finally, for a given induced dipole d c the ground-state of an ensemble of polar 
molecules is a crystal for mean interparticle distances a m i n < a < a m ax, where 
a max = d 2 m/h 2 r c corresponds to the distance at which the crystal melts into a 
superfluid. For SrO (RbCs) molecules with the permanent dipole moment d c = 8.9D 
(d c = 1.25D), 

a min ~ 200nm(100nm) , while a m ax can be several fim. Since for 
large enough interactions the melting temperature T m can be of order of several /iK, 
the self-assembled crystalline phase should be accessible for reasonable experimental 
parameters using cold polar molecules. 



4. Specific implementations with atoms and molecules in dipolar crystals 

In this section we consider a mixture of two species of particles confined to one 
dimension. The first species of particles comprises (strongly interacting) molecules 
with dipole moment d c , forming a one-dimensional crystal. The second species of 
particles can be either atoms [see figure Hfb)] or molecules of a second species, with 
dipole moment d p -C d c [see figure [He)]. The former interact with the crystal 
molecules via a short range pseudopotential proportional to an elastic scattering length 
a cp , while the latter interact with the crystal molecules via long-range dipole-dipolc 
interactions. For both configurations, we obtain explicit expressions for all parameters 
characterizing the coherent and the incoherent dynamics in the system. These one 
dimensional setups can be readily generalized to two dimensions. 
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Figure 3. We show the dispersion for a ID and 2D dipolar crystal as a function 
of the quasimomentum in units of the crystal recoil energy -Er, c - In ID (a) the 
dispersion is peaked at the zone border and tends to zero oc q for small momenta. 
The 2D dispersion has two acoustic branches, a longitudinal and a transversal 
one. We plot the two dispersions as a function of the quasimomentum, choosing 
a path in the first Brillouin zone that is outlined in the inset of the figure. 



4--1- Neutral atoms moving inside a crystal tube 

As a first realization, we consider a setup, where an ensemble of neutral atoms is 
confined by an optical trap to the same ID tube as the dipolar crystal, see FigurefTJb), 
say along x. For simplicity we assume the trap for the neutral atoms and the polar 
molecules to have the same harmonic oscillator frequency u±. 

An atom and a molecule inside the tube interact via a short range potential, 
which we model in the form of an effective ID zero range potential, 

V cp (x - X) = g cp 5(x - X) 

Assuming that the 3D scattering length a cp is (much) smaller than the harmonic 
oscillator length of the traps, a cp <C a p ,j_ = (h/m-pLn^) 1 / 2 , the effective ID coupling 
strength is given by g cp ~ 2hiu±a cp . In the following we focus on positive scattering 
lengths, a cp > 0, corresponding to a situation where the atoms and molecules 
effectively repel each other, cf. g cp > 0. 

4- 1.1. Tight binding limit and Hubbard models. For a "frozen" crystal, i.e. without 
phonons, the molecules are at their cquilibirum positions, ja + a/2, which provide for 
a static periodic potential for the atoms, 

V P {x)=g cp Y i S(x-ja-a/2). (41) 

3 

The dynamics for a single neutral atom is then determined from the static 
Hamiltonian H p — p 2 /2m p + V p (x), corresponding to the Kronig-Penney model with 
a potential comb of strength g cp . Its energy spectrum is given in the form of a band- 
structure, E n ^ q (with band-index n = 0, 1, . . .), which is obtained from 

JLy cot ( m± ^\ = ^± W ith a^i^A 1 ' 2 (42) 
4" ^ \ 2 ) g cp \E R , P J 1 ' 
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(a) u/e Rj ,a/e Rj , (b) Mq [NVVJ 4 K gce m 




Figure 4. (a) The bandwidth 4J of the lowest band (solid blue line) and the 
gap A to the first excited band (solid red line) for a neutral atom scattered from 
a ID potential comb of strength g cp and lattice spacing a. All energies are given 
in terms of the particle recoil energy i?R lP - The dashed line denotes the coupling 
strength, g cp /a t*£ Er iP /2, where the gap and band-width are equal, (b) The 

corresponding particle phonon coupling M q in units of rV^-s/NE^p as a function 
of the quasi momentum q of the atom. The interaction is linear and peaked for 
large q, while it shows a square-root behavior for small quasi momenta (see text). 

where E^ tP = H 2 tt 2 /2m p a 2 denotes the recoil energy of an atom. We denote the band- 
width of the lowest band by 4 J = Eq^u — Eo t o = -Ert.p — E n / a , and the gap to the 
first band by A = Ei t „/ a — Eq^/o. — Ei,ir/a — Eh,p- The latter are both shown in 
figure[4|a) as a function of the coupling strength g cp . We notice that for g cp > E^ p a/2 
(indicated by a vertical dashed line), the gap exceeds the band- width, A > 4 J. In 
the tight binding limit, cf. g cp ^> E^ p a, the dispersion relation in the lowest band 
becomes E 0iq sa 4Jsin 2 (ga/2), with J ps 2E^ p /n 2 g cp + 0(aE RiP /g cp ) and A rj 3S RiP . 
The Wannier-functions for a particle become localized at site j and are approximated 
by Wj(x) i=s cos[7r(a; — ja)/a]/ \fa~j2 for \x — ja\ < a/2 and zero otherwise. 

Obtaining the tight-binding limit requires that the ratio, g cp / 'aE^ p = 
2a cp a/ir 2 a 2 ± (largely) exceed the value « 1/2. We notice that for current state-of the 
art optical traps, one can achieve harmonic oscillator lengths as small as a p .± ~ 20 nm, 
and taking a "typical" 3D scattering length of a cp ~ lOOao ~ 5 nm, we get that 
g C p/aE^ p > 1/2 is attained for lattice spacings a > n 2 a 2 ± /a cp ~ 200nm. 

Analogous to the atom-molecule interactions, we model the interactions between 
two neutral atoms by a contact potential with a coupling strength given by (? pp 
h\j_\_a pp for the 3D atom-atom scattering length a pp <^ a±. Then in the tight-binding 
limit the atom-atom interactions reduce to repulsive onsite energy shifts only 

Vij = g pp J dxwi{x) 2 Wj{x) 2 sa ^^^j- 

The dynamics for an ensemble of bosonic atoms in the crystal is then described by 
a single band Bose-Hubbard model with hopping rate J and onsite repulsion Vu, 
provided that Vu <C A. On the other hand for an ensemble of (spin-polarized) 



Quantum Simulations of Extended Hubbard Models with Dipolar Crystals 



16 



fermionic atoms, we notice that the system reduces to a lattice model with hopping 
rate J and no interactions. 

The atoms couple to the crystal molecules via a density-displacement 
interaction jH^|. To first order in the displacement (see | Appendix B.2 |, the coupling 
constant reads 

M g = ( h Wcp(g) = —\{w^77W« ( 43 ) 

\2N c m c LU q J a y Nm c uj q 

where V cp is the Fourier transform of the atom-molecule interaction potential, and f3 q 
is the Fourier transform of the square of the Wannier-functions. For g cp ^> aE^ p the 
latter is 

(3 q = [ dxwaixfe** * *? ^ f 
J A-K z qa — q 6 a 6 

The latter decreases with increasing q from (3q = 1 to /3 n / a — 8 /Sir w 

0. 85. The (monotonical) dependence of the coupling constant M g on the quasi- 
momentum is shown in figure H[b). In particular, it has a maximum, ~ 
(8g cp /3a 2 )(2h/Nm c LJY)) 1 / 2 , at the band-edges, while for small quasi-momenta q it 
vanishes as IqI 1 / 2 , i.e. M q w (2 3cp /a 2 )(n|qa|/7VTO c w D ) 1 / 2 + 0(ga) 5 / 2 . 

Finally, let us address the validity of the single band approximation in the 
Hubbard model {I), when coupled to phonons. For simplicity let us consider the 
limit of vanishing interactions Vu — and a weak coupling M q : We notice that, 
for Hujd < 4J + A the second band is gapped from the branch of acoustic phonons, 
and therefore higher band excitations are (strongly) suppressed and can be neglected. 
Since -Er, p < 4 J + A < 31?r iP , this requires a mass ratio m p /m c < 3/^2rd- While 
this for a soft crystal with rd ~ 1 merely implies that m p < m c , for a stiff crystal with 
rd ~ 200 this requires to p < 0.15m c , which is quite restrictive. In the latter case, that 
is for a stiff crystal and comparable masses, cf. ^/2rdm p /3m c > 1, we notice that the 
first excited band would cut the phonon branch at a frequency ~ 4J + A. However, 
in this regime a single band model may still hold, if one restricts the initial phonon 
population to sufficiently low temperatures, i.e. for k^T <C hA. In the tight binding 
limit this requires temperatures T <C (3m c / ' \f2rdm p ) x hwa/k^ which even for a stiff 
crystal with rd ~ 200 and comparable mass ratio m p ~ m c yields temperatures on the 
order of Q.lhuj^/k^,. These are smaller than the melting temperature of the crystal, 
T m ~ 0.1y/2rdh~oJB /k-g, (in ID), and are thus reasonable, even for a "soft" crystal with 
r d — 1. 

4..1.2. Extended Hubbard model for atomic polarons inside a dipolar crystal As we 
have seen in Section 12.11 it is convenient to change from a picture of bare atoms and 
crystal phonons, to one of atoms dressed by their surrounding crystal displacements, 

1. e. polarons. The corresponding displacement amplitudes u q for the dressing then 
are 

u M± J 2 * 2JJ » 3 -22- (44) 

U " to, VF(186C(5))3/^« w 3/ 2 ' (M) 

which at the band-edge attain their minimum, u^/ a « imU^/N 1 ' 2 , while they diverge 
at small quasimomenta as ~ 4.94|ga| -1 / 2 . In (|44p we introduced the coupling ratio 

Sep 5c P m c 1 



U 
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(b) V u /E RtP ,J/E Kj , 
0.2- 
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Figure 5. (a) The ratio S T /U$ as a function of the (dimcnsionless) temperature 
k^T/hwY) with the dimensionless prefactor Uo = (g cp / aE^^ p ) / (m-p / m c )r^J A . The 
dashed line denotes the value St=o/Uq and we find that for small temperatures 
kftT/hu>Y) <C 1 the influence of the temperature on St is negligible. St depends 
quadratically on Uq which realistically takes up values inbetween 10 -2 and 10. 
Thus we switch between the weak (St <S 1) and strong (St 3> 1) coupling 
regimes by choice of the dimensionless prefactor Uo- (b) The full particle- 
particle interaction Vij between two extra particles at intcrparticle distances 
I* — j\ = D) lj 2, 3 for r d = 100 and m p /m c = 0.1 as a function of the coupling 
g cp in units of the particle recoil energy. Notice that the sign of the interaction 
alternates with every site and that the total interaction strength decreases with 
the interparticle distance as oc 1 / 1 % — j \ 2 . 



which increases linearly with the atom-molecule coupling constant g cp /aE^ p and the 
mass ratio of m c /m p , but is inversely proportional to the "stiffness" of the crystal rd- 
The dependence of Uq x m p /m c = g^/o-Ecpf^/ 4 on the coupling g cp /aE cp and r c i is 
also shown as dashed contour lines in figure [B^a). We see that for, e.g., a mass ratio 
of m p /m c ~ 1 (e.g. for a gas of Cs atoms in a crystal of LiCs molecules ) Uq can 
take values ranging from Uq rs 0.01 (at g cp /aE^ p = 1/2, rd = 200) to Uq w 5 (at 
gcp/aEft p = 5, = 1), while still having a crystal and being in the tight binding 
limit. For a mixture with a mass ratio of m p /m c ~ 1/20 the coupling Uq is significantly 
larger, i.e. Uq » 50 (at g cv /aE-s.. v = 5, rd = 1 for m p /m c ~ 1/20). 

Due to the dragging of the surrounding phonon cloud, the hopping rate for the 
polarons, J, compared to the (bare) hopping rate, J, is suppressed by the factor 
where the exponent is, cf. (fTl) . 



J/J 



-S T 



St — 



8tt 4 U§ 



(186C(5)) 3 / 2 N Y w; 



7-3 ^ sin ' 



21, b t) ' 



(45) 



We notice that St increases quadratically with Uq (cf. the coupling <? C p), whereas the 
ratio St/Uq depends only on the temperature T (in units of the Debye frequency), 
which is shown in figure Efa). In particular we remark that for T = the exponent 
scales with the coupling ratio Uo as St=q ~ 0.9C/q [indicated by a horizontal dashed 
line in figure |S[a)] , while St/Uq only weakly increases with the temperature T. 
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Figure 6. (a) Contour plots of the coupling ratio times the mass ratio 
Uq X (m p /m c ) (dashed contour lines) and of the ratio of the (bare) atomic 
tunneling rate and the Debye frequency of the crystal J/tko-Q (solid contour 
lines) as a function of the atom-molecule coupling g cp / 'ciEr p and the stiffness 
of the crystal r^. The ratio Uo(m p /m c ) determines the overall strength of the 
displacement amplitudes u q for a polaron, while the ratio J/tjuJu characterizes the 
serparation of crystal and interaction time (see text), and serves as "the" smallness 
parameter in the derivation of a master-equation (1131 , (b) The spectral density 
for the swapping of two particles on neighboring sites, Jq^(lu), as a function of 
the frequency lu/ujd. It displays a Van Hove singularity at the Debye frequency 
u>d where it diverges as ~ (u>rj — a;) -1 / 2 . 



The phonon-coupling provides phonon mediated particle particle interactions of 
strength 

y(i) = _ 4tt 2 9% (m c /m p ) ^ £_p , aaH _ j]) m 
» 93C(5)£ R . P Nr d ^ W 2^ C0S W 0\h W 

for two polarons at sites i and j, respectively. We notice that are temperature- 
independent and vary in sign and magnitude with the separation i — j, i.e. they are 
attractive for even i — j and repulsive for odd i — j, while their absolute value decreases 
with increasing inter-polaron separation \i — j\. In figure [H^b) we show the leading 
contributions for the total off-site shifts, which are purely induced by the phonons, 
Vij = V^j for i 7^ j, as a function of g cp /aE^ p , cf. the leading off-site terms decay 
as Vij/2E P w 0.16/|* — j\ 2 . In addition we also plot the corresponding (modified) 
hopping rate J for zero temperature T = 0. We notice that near g c -plaE^ p ~ 1.6 
the nearest neighbor interactions become comparable with the effective hopping rate, 
y%,i+x ~ J- 

For bosonic particles, the phonon-mediated interactions also include an attractve 
on-site shift, the value of which is exactly twice the polaron shift, 

V£ ] = -2E P w -1.54 x U^huj u . 

This can (in principle) lead to a collapse of the system, as it favors the piling up of 
polarons at a single site. However, since the full interactions comprise the bare and 
the phonon-mediated interaction, the total onsite shift, Vu — Vu — 2E p , is positive for 
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Va/2 > E p , which we require to ensure the stability of the bosonic system. We remark 
that by resorting to tune g pp via a Feshbach resonance, one can tune the onsite-shift 
up to Vu ~ A, without breaking the single-band approximation in the Hubbard model 
fl}. Thus we notice that in principle the stability of the system can be guaranteed, 
provided E p < (Vu/2) < A/2. 

4- 1.3. Corrections to the extended Hubbard model In the following we are interested 
in higher-order corrections to the effective Hubbard model H$ of © , which we derived 
in Section [2] in terms of the spectral densities J^J (u>) for correlated nearest-neighbor 
hopping events i — ► j and k — ► I. For our atomic-crystalline mixture we find from (|20[) 
that the latter are given by 

J « {W) (186C(5))3/2 ( w / WD )3 A /^3^2^ ^ ^ 7 > 

where we took (3 q ps 1 and oj q ps wd sin(ga/2) and q u = aicsiniuj / loy>) / a. The spectral 
density for the "swapping" of two particles on neighboring sites, Jqi(u>) is shown in 
Figure BJc) and shows a Van Hove singularity at ui — ► Ud, due to the l/y/uiu — oj 
divergence of the density of states for the crystal phonons. 

Strong coupling limit St 3> 1- The value St is determined from equation (|45[) . 
Values St 3> 1 are obtained for large coupling ratios Uq and/or high temperatures 
T. However, already for a (reasonably small) ratio Uq > 1.1 we have St > 1 
at T = and thus we are in the strong-coupling regime for all temperatures 
T > 0. In this limit the main corrections to the Hubbard model © are due to 
"swap" processes. The corresponding rates and coefficients for St 3> 1 are well 
approximated by the expressions (|226[) - (|22<i[) . which are shown in figure [7] and figure [8] 
as a function of the coupling ratio Uq and the temperature T. Notice that in (|22oj) - 
([22d|) the parameter B = 4tt 6 [/ 2 /3(186C(5)) 3 / 2 ps 0.48J7 2 while the ratio A T /B exceeds 
(k B T/huj D ) tanh(fiw D /2fc B T). 

In figure [T^a-d) we show the leading contributions in the strong-coupling regime 
as a function of the ratio Uq and the temperature T. In particular, panels (a), (b), (c), 
and (d) are contour plots of the incoherent rate HTqI/J, the coherent shift Aq1/J, 7 
and 5, respectively, as obtained from the strong coupling approximation (|22aj) - (|22rij) . 
The two dashed lines in each panel signal where St = 10 and St = 1 ; and the strong- 
coupling approximation is valid (St > 1)- We remark that in panels (a,b) T/J, A/ J 
are divided by the (small) ratio J/Hlod, while in (c,d) 7 and 5 are divided by the even 
smaller ratio (J/Hcod) 2 . We notice that in the spirit of the master equation approach, 
all quantities in the figures are plotted at finite temperature. The figure shows that, 
in a wide range of parameters, corrections to the coherent time evolution determined 
by Hs can be made small in the strong coupling regime. 

Weak coupling limit St <C 1 ■' In Figure [8] we show the leading corrections to 
the extended Hubbard model Hs as a function of the coupling ratio Uq and the 
temperature T, as obtained from the weak-coupling approximations (|27aj) - (|27d|) - In 
particular, the solid lines now indicate contours of the largest value attained for (a) 
T/J, (b) A/ J, (c) 7 and (d) S. In panels (a,b), T/J, and A/ J are divided by the small 
ratio J/huj-Q, while 7 and 5 in (c,d) are divided by the even smaller ratio (J/Swd) 2 - 
The two dashed lines indicate St — 0.1 and St = 1, where the latter delimits the 
range of validity of the weak-coupling expressions (|27a|) - (|27d|) (e.g., Uq < 1 for T = 
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Figure 7. Leading corrections to the extended Hubbard model for a atomic- 
crystalline mixture in the strong coupling regime. St 2> 1, as function of the 
temperature of the crystal T and the coupling ratio Uq: The solid lines indicate 
contours for (a) the incoherent "rate" rJJ and (b) the coherent "shift" AJJ (in 
units the effective hopping rate J), (c) the incoherent coefficient and (d) the 
coherent coefficient <5qJ as obtained from the strong-coupling approximation 1 122 al - 
l!22d . respectively. The two dashed lines in each panel represent the contours 
where St = 1 and St = 10, respectively. They designate the area, where St > 1, 
and thus the strong coupling approximations holds. Notice that in (a,b) T/J, 
A/J are divided by the small ratio J/Hui-q, while in (c,d) 7, 8/J are divided by 
the even smaller ratio (J/Tjud) 2 - 



while Uq < 1/2 for T = hu-o). We notice that in the area where St ;$ 0.1 the ratios 
shown in panels (a-d) are smaller than « 1, and thus all corrections are strongly 
suppressed compared to J, provided J <C hu>. 



4-2. Polar molecules interacting with a one- dimensional crystal 

As a second configuration, we consider a setup where polar molecules of a second 
species are trapped at a distance b from the crystal tube, under one-dimensional 
trapping conditions [see figure H{c)]. An external electric field aligns all dipoles in the 
direction perpendicular to the plane containing the two tubes. Molecules trapped in 
the two different tubes interact via long-range dipole-dipole interactions. 



4-2.1. Tight binding limit and Hubbard models For crystal molecules fixed at the 
equilibrium positions with lattice spacing a, the particles (that is, the molecules of the 
second species) feel the following periodic potential 

V cp {x) = ^rJ2 [(6/a) 2 + (x/a-i-l/2) 2 ]3/ 2 ' (48) 
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Figure 8. Leading corrections to the extended Hubbard model for an atomic- 
crystalline mixture in the weak coupling regime, St C 1, as function of the 
temperature of the crystal T and the coupling ratio Uo: The solid lines indicate 
contours for the largest (a) incoherent rate r max and (b) coherent shift A max (in 
units the effective hopping rate J), (c) incoherent coefficient 7 max and (d) coherent 
coefficient <5 max as obtained from the weak-coupling approximation H27aM27crt . 
respectively. The two dashed lines in each panel represent the contours where 
St = 0.1 and St = 1, designating the area, where St < 1, and thus the weak 
coupling approximations hold. Notice that in (a,b) T/J, A/ J are divided by the 
small ratio J/hui-o, while in (c,d) 7, S/J are divided by the even smaller ratio 
(J/^d) 2 . 



where d p is the induced dipole moment of the second-species molecules. The potential 
above has a depth 

V = V cp (a/2) - V cp (0) ~ v e~ 3b / a E R , p /(b/a) 3 , 

which determines the band-structure for the particles, with 

va = {d p /d c )(m c /m p )r d , (49) 

and i?R,p = h 2 ir 2 /2m p a 2 the particle recoil energy. The lattice depth Vq is 
shown in Fig. figure O^a) to have a comb- like structure for b/a < 1/4, since the 
particles resolve the individual molecules forming the crystal, while it is sinusoidal 
for b/a > 1/4. Figure \§0.b) shows the width 4 J of the lowest-energy band, with 
J/E R , P ~ {Vq /£ R ,p) 3 / 4 e- V v 'o/-Br,p f or / a > together with the energy gap 

A ~ (4Vbi?R,p) 1 ^ 2 , as a function of b/a and for vq = 1 and 50. For a single particle, 
the single-band model is valid for 4 J < A. Figure [He) is a contour plot of the regimes 
of validity of the single-band model as a function of b/a and vq- 

When more particles are considered, the strong dipole-dipole repulsion between 
the particles acts as an effective hard-core constraint [37]. We find that for 4 J < A 
and dp <C d c the bare off-site interactions satisfy ~ d p /(a\i — j |) 3 < A, and thus 
the single-band model is still valid. 
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Figure 9. In panel (a) we plot the potential depth Vb as a function of b/a for 
different values of in units of the crystal recoil energy to give an idea about the 
extra-particle crystal interaction potential which determines the bandstructure. 
The hopping amplitude 4J and the gap to the first excited band A are shown 
in panel (b) as a function of b/a in units of the particle recoil energy. The 
bandstructure strongly depends on the ratio of the particle-crystal interaction 
energy over the kinetic particle energy vo = (d p / d c )r,i/ {m p /m c ). The hopping 
amplitude decreases while the gap increases rapidly with increasing coupling 
strength and vq. The single band model is only valid where the gap exceeds 
the bandwidth A > 4J. 



In this configuration, the particlc-phonon coupling as obtained from equa- 
tion (|A.13j) is given by 

where /Ci denotes the modified Bessel function of the second kind, and f3 q = 
J dxe lqx \wo(x)\ 2 , with wq(x) the lowest-band Wannier functions. Figure [TUJa) shows 
that for b/a small-enough, such that the single-band approximation is fulfilled for all 
■Do [see figure [9jc) above], M q becomes peaked at large quasimomenta q. We notice 
that consistency with the requirement of a stable crystal implies that the variance of 
the fluctuations of the crystal molecules around their equilibrium positions induced 
by the presence of a particle localized at a site j, (Svij), be small compared to the 
interparticle distance a, that is (5vij)/a < 1. For a given ratio d p /d c , this limits 
how small the ratio b/a can realistically be, in order to avoid that the inter-species 
interactions destroy the crystalline structure [571 ] . For example, for a ratio d p /d c w 0.1 
the ratio b/a can be as small as b/a ~ 0.2. 



4-2.2. Extended Hubbard model for molecular polarons inside a dipolar crystal 
We continue by determining the modified Hubbard parameters J and Vij for this 
configuration. Here, the parameter St, which determines the regime of interactions, 
is given by 

s = Z2r 1 J 2 {d v /d c f 1 M 4 /C?(%|) a2 ^ 2 fqa^ _ u ( hw q ^ 



(186C(5)) 3 / 2 (b/a)iw 3 q q V 2 V 1 \2k B T) K ' 
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Figure 10. Figure (a) shows the dependence of particle-phonon coupling M q 
on the quasi momentum qa for different values of b/a. Notice that as expected 
the coupling strength increases with small b/a and large vo- For small q it tends 
to zero like q 1 / 2 /(b/a) 2 . Figure (b) shows how the spectral density behaves as a 
function of w for different values of b/a on a logarithmic scale. It depends strongly 
on the ratio b/a and tends to zero like ui for small frequncics. 



The latter depends strongly on the ratio b/a and is proportional to r^ 2 (d p / d c ) 2 . For 
a given and d p /d c ratio, the regimes of weak and strong coupling, St *C 1 and 
St ^ 1) respectively, can be directly determined from figure [TT1 which is a contour 
plot of St as a function of the dimensionless temperature k^T/hw^, and the ratio 
b/a. As in the previous model St increases with increasing temperature and particle- 
phonon coupling (that is, with decreasing ratio b/a). 

The phonon mediated interaction as determined from equation JSJ is given by 



r(1) _ 16r d (d P K) 2 1 v (qa^qiM) 

93C(5)tt 2 R ' c N^ {b/a) 2 w 2 Mq 



vx' - i: (52) 

As in the previous model, the phonon mediated interactions show oscillations which 
for b/a < 1/4 decay slowly as l/\i — j\ and are thus long-ranged. Depending on 
their sign, they can enhance or reduce the bare dipole-dipole repulsions between two 
second-species molecules. The phonon-mediated interaction is strong and dominates 
the full particle-particle interaction Vij for small b/a. This is shown in figure [TTT bl 
which is a plot of Vij /Vij as a function of b/a, where for b/a < 0.4 the value of Vij /Vij 
can even change sign. With increasing intertube distance, the phonon-mediated term 
becomes small compared to the bare value Vij . 

4-2.3. Corrections to the extended Hubbard model In the following we arc interested in 
coherent and incoherent corrections to the time evolution determined by the effective 
Hubbard Hamiltonian iJg of ([3]). These can be calculated in terms of the spectral 
density, which has the form 

T M M 64(d p /4) 2 ry 2 0,3 MKijbq^Pl 

J ^ [W) - 7r(186C(5))3/2 W 3 (b/aY^l^ 9 ^ (qwh {b6) 
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Figure 11. The quantity St depends on the temperature k^T, the ratio b/a and 
is proportional to r^J 2 (d p /d c ) 2 . In figure (a) we show how St behaves with the 
dimensionless temperature k^T/hui^s and the ratio b/a. By fixing the dipole ratio 
and r d we can determine where the weak and strong coupling regimes are valid. 
Figure (b) shows the full particle-particle interaction Vij, the sum of the bare 
dipole-dipole repulsion with the phonon-mediated interaction, in units of the bare 
nearest neighbour interaction Vj t+1- The full interaction decays with increasing 
distance and shows an alternating sign for small ratios b/a where the phonon- 
mediated interaction J dominates. In figure (c) we show a contour of the 
separation of bath and interaction timescales, J/ftujD, that play an importent role 
when discussing the validity of our model and the amplitude of the perturbative 
corrections. 



where we have used uj q » ujb sin(ga/2), and q w a — 2 arcsin(u> / Wd) . The spectral den- 
sity tends to zero linearly for small w and shows an integrable (van Hove) singularity 
at w = lub- The spectral density Jqi(w) is plotted in figure [TOl for a few values of b/a. 

General expressions for the corrections A g (T), T q (T), 7 g (T) and S q (T) are given 
in (|17a[) - (|17rfp . For the configuration that we consider here, they depend on b/a, the 
temperature ks,T and the dimensionless parameter (d p /d c ) 2 r^/ 2 . The quantities T q (T) 
and Aq(T) are proportional to the ratio J/huii), while S q (T) and Jq(T) are propor- 
tional to (J/Hujb) 2 - Thus, for later convenience, in figure [TTTc) we plot J/hwu as a 
function of b/a and r^ for a realistic choice of the dipole and mass ratios, d p /d c = 0.2 
and m p /m c = 0.5, respectively. The figure shows that for this choice of parameters 
the ratio J/Swd is (much) smaller than one for all plotted values of b/a and ra, and 
in particular it is e.g. of order ~ 10 -4 for reasonable values b/a w 0.4 and w 150. 

Strong coupling limit St S> 1: Here we are interested in giving examples of 
the importance of the corrections for realistic parameter regimes, compared to the 
characteristic energy J of the polaronic Hamiltonian Hs- Thus, in figure fT2l we show 
contour plots of the quantities nrJ°(T)/J, Aj?(T)/J, 7^(T) and 6${T) as a function 
of b/a and the dimensionless temperature k^T/tuu^, and the ratio {d-p/dc) 2 ^ 2 . 

As explained in Section 2, these quantities correspond to the corrections for 
the "swap" process (i — I and k — j), which is not suppressed exponentially by a 
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Figure 12. Ratios of the only non-negligible corrections that contribute to the 
r o (T) » 2T™(T), A,(T) « 2AJ?(T), 7 ,(T) « 137$ (T) and 



S q (T) -C 12<5qJ(T) over the effective tunneling rate J in the strong coupling 
limit for a single extra particle. These ratios are shown as functions of b/a, the 
dimensionless temperature k^T/huiY) and (dp/d c ) 2 r^ 2 in two distinct plots where 

we first fix (d-p / d c ) 2 r X J 2 = 0.2 and then the temperature as iflT/fiUD = 0.1. The 
blue dotted lines show the value of St and the strong coupling approximation 
breaks down below St = 0.5. When including the separation of timescales J/tkJD 
that is shown in figure lllf c) we find that the corrections ■jq (T) and & q (T) are 



negligible values for large St and high temperatures. 



an d is therefore the dominant correction in the strong-coupling 
In particular, we can estimate A q (T) s» 2Aj"(T) and 



factor oc exp( 
limit St 3> 

T q (T) w 2T^((T), while upper bounds for 6 q (T) and j q (T) can be estimated as 

u max 

(T) < 125$ (T) and (T) < 12j$(T). 
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Figure 13. Ratios of the maximal eigenvalues of the corrections to the master 
equation over the effective tunneling rate J in the weak coupling limit, St ■C 1, 
for a single extra particle. They plotted as functions of the ratio b/a, the 
dimensionless temperature ksT/hiLij^ and (d p /d c ) 2 r^ in two distinct plots where 

we first fix (dp/dc) 2 ^ 2 = 0.2 and then the temperature as k^T/tiWY) = 0.1. The 
blue dotted lines outline the value of St and the wek coupling approximation 
breakes down for values of St > 0.5. When taking the separation of timescales 
into account that is outlined in figure 1111 we find that for St < 0.5 all the 
corrections to the master equation are negligiable compared to J. 

Panels (al), (bl), (cl) and (dl) of figureQJshow results for hT$(T)/ J, Aj?(T)/J, 
7qi(T) and 5q® (T) as a function of k^T/huj^,, respectively, while the ratio (d p /d c ) 2 r^ 2 
is fixed to the reasonable value 0.2. Panels (a2), (b2), (c2) and (d2) show the 
corrections as a function of (dp/d c ) 2 r^/ 2 , with the temperature fixed to the value 
k B T/huj D = 0.1. 
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Figure [T2Ta.lV(d. 2) show that, for reasonably small J/Swd [see figure UTTc)]. 7g° 
and tend to remain small in the range of shown parameters at finite T. However, 
the rate hTl° [panels (a.l)-(a.2)] and the energy [panels (b.l)-(b.2)] can exceed 
the effective tunneling rate J when the latter is strongly suppressed for strong cou- 
plings > 1. While large fiTj-jJ can in principle lead to significant decoherence, 
and thus a transition from coherent hopping to thermally-activated hopping for large 
enough temperatures, we find that for reasonable temperatures hsT/hun < 0.1 these 
processes are strongly suppressed. On the other hand, the self-energies Aq°(T) can 
lead to significant modifications to the coherent-time evolution determined by Hg by 
providing next-nearest neighbor hopping and sizeable off-site interactions in the strong 
coupling regime St 3> 1. 

Weak coupling limit St *C 1: Figure 1131 show the corrections to the coherent- 
time evolution given by H§ as a function of b/a and the dimensionless tempera- 
ture kBT/hwD [panels (a.l),(b.l),(c.l) and (d.l)] and the ratio (dp/ 'd c ) 2 r^ 2 [panels 
(a.2),(b.2),(c.2) and (d.2)], in the regime of parameters where the single-band approx- 
imation is valid. We find that the weak coupling limit covers most of the accessible 
parameter regime for b/a. Figure [T3] shows the maximal eigenvalues fir max , A max , 
7max and <5 max for all four corrections. The central result here is that we find that 
the latter are negligible compared to the coherent hopping J in the region where the 
weak coupling expansion is valid. In the figure, as a reference to identify how good 
the "weak-coupling" expansion is we also plot the values of St for the various regimes 
of parameters. 

Finally, we conclude this section by providing an example of the regime of validity 
of our model configuration. We consider a crystal of SrO where second-species 
molecules are KRb. The dipole and mass ratios are d p /d c w 0.08 and m p /m c w 1.2, 
respectively. We find that our treatment properly accounts for the system dynamics 
for separations 0.2 < b/a < 0.7, provided > 80. 

5. Conclusion 

In this work we studied the realization of lattice models in mixtures of cold atoms and 
polar molecules, where a first molecular species is in a crystalline configuration and 
provides a periodic trapping potential for the second (atomic or molecular) species. 
We have treated the system dynamics in a master equation formalism in the Brownian 
motion limit for slow, massive, particles embedded in the molecular crystal with fast 
phonons. In a wide regime of parameters the reduced system dynamics corresponds 
to coherent evolution for particles dressed by lattice phonons, which is well described 
by extended Hubbard models. For two realistic one-dimensional setups with atoms 
and molecules these lattice models display phonon-mediated interactions which are 
strong and long-ranged (decaying as l/\i — j\ ). The sign of interactions can vary 
with distance from repulsive to attractive, which can possibly lead to the realization 
of interesting phases of interacting polarons in one dimensional dipolar crystals. This 
study, and extensions to two dimensions will be the subject of future work. 
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Appendix A. Hamiltonian for particles moving in a crystal 

It is the aim of this section to derive the Hamiltonian |T]) starting from a generic 
mixture of two interacting species of atoms or molecules. 

Appendix A.l. Effective continuum Hamiltonian 

A mixture of two interacting species is confined to one or two dimensions by a strong 
optical trapping potential. In [52| it is shown how such a trapped system can be 
reduced to an effective lower dimensional model in the low energy limit. The effective 
one- or two-dimensional Hamiltonian reads 

H cB = H c + H p + H cp . (A.l) 

Here H c describes the effective motion of the crystal particles, H p the dynamics of the 
extra particles with an extra particle- crystal particle interaction given by H cp . From 
Hamiltonian (|A.1[) we identify 

H cp = Y,\ ^r, R, ;. 

i,3 

where we denote pj(r<) and P,;(R,) the momentum (position) of extra particles and 
crystal particles with masses m p and m c , respectively. The sums range over all the 
respective particles in the mixture. The interaction between two particles from the 
species a, [3 € {c, p} in a distance r from each other is denoted by V a p(r). 

We take the continuum Hamiltonian (|A.1|) as the starting point for the following 
discussion and derive our lattice model from it. 

Appendix A.2. Crystal Hamiltonian 

In the crystalline phase the particles are characterized by small fluctuations of their 
positions around their equilibrium positions R^. Thus we expand the potential in a 
Taylor series to second order in the displacements Uj = R 3 — R° about the equilibrium 
positions [H^ as 

V cc (Ri - Rj) w U CC (R° - R°) + u.D. u ,. (A.5) 

with the tensor 

D tj = -V ® VV CC (R° - R°), (A.6) 

which is readily diagonalized and provides the dispersion relation u)q,\ for the phonons 
with quasimomentum q and polarization e^. We express the displacement in 



- Rj), (A.2) 

- r,), (A.3) 

(A.4) 
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terms of the bosonic creation (annihilation) operators a q A ( a q,\) f° r the corresponding 
phonons as 



q.A V M ' 

where N denotes the total number of crystal particles. We write the Hamiltonian as 
H c = 2J ^q,Aa qi A a q,A (A. 8) 

q,A 

and can identify m c Wq A from the eigenvalues of Dij from which we are able to read 
off the phononic dispersion relation, see | Appendix B| 

Appendix A. 3. Extra particles inside the crystal 

The dynamics of extra particles inside the crystal and their interaction with the crystal 
is described by H p + H cp . 

In analogy to the proceedings of the last section we are interested in the dynamics 
of a stiff crystal and therefore we expand the particle-crystal potential in the small 
displacements of crystal constituents about their equilibrium positions, as 

V p (r) = ]T V cp (r - Ri) » £ V cp (r - R°) + ]T u,Vy cp (r - R°). (A.9) 

i i i 

Here we have introduced the potential V p (r) that an extra particle at position r feels 
from the entire crystal. The lowest order provides a static periodic potential 

^ P (0) (r)-E^p( r - R ")- 

i 

We note that in writing equation (|A.9p wc only retain the first non- vanishing correction 
to the static trapping potential, i.e. the term linear in the displacement u,. Thereby 
we neglect terms of second (and higher) order in Uj, which are expected to (merely) 
provide a renormalization of the phonon-spectrum, i.e. when including those terms in 
equation (|A.6|) . 

The extra particles that we consider are confined to a plane/tube parallel to a 
crystal plane/tube, as pictured in figure[TJ Hamiltonian (|A.3|) together with the zeroth 
order contribution to the particle-crystal interaction, 

ff^ifp+XVWfc), (A.IO) 

i 

describes interacting particles for which the entire crystal provides a static periodic 
trapping potential, the ingredients for a simple Hubbard model. The particle crystal 
interaction then determines the band structure. 

We consider a single band model where extra particles cannot be excited to the 
second band that is separated from the lowest band by an energy gap A. Then in 
the low energy limit the extra particles localize at single sites of position r° and their 
wave functions become Wannier functions of the lowest band u>o(r — r,). In second 
quantization the field operator can then be expanded as 

^p(r) =^c 4 w (r-r°), 
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where c\ and Ci denote the creation and annihilation operators of extra particles at 
site i, which obey the canonical bosonic (fermionic) commutation (anticommutation) 
relations for bosons (fermions) 



For such a model, see [76[, the dynamics of the extra particles is described by 
the hopping amplitude J that calculates from the overlap of the wavefunctions at two 
neighbouring sites and the interaction between two particles located at sites i and j 
is given by 

Vu « J dvdv'\w {v - r°)| V pp (r - r')K(r' - r°)| 2 . 
Hamiltonian (|A.10[) can then be written as 

H ' P = - 3 E ^ + \ E Vij&h *- ( A - n ) 

The single band approximation is valid if all particle energies are (much) smaller than 
the gap, cf. J, Vij <C A. We remark, that the Debye frequency fouin in our models is 
typically (much) larger than the gap while the particle phonon coupling is dominated 
at high frequencies, see discussion below. In order to avoid excitations beyond the 
gap by the coupling we put a constraint on the temperature in a way that all phonon 
modes with energies larger than A, are essentially unoccupied. 

Let us now focus on the higher order terms of the interaction in equation (|A.9l) 
which correspond to the backaction of the crystal on the extra particles. The remaining 
part of the Hamiltonian is given by 

i 

and describes a dynamic coupling of the particles to the vibrations of the crystal. In 
second quantization for the extra particles the remaining Hamiltonian is obtained (cf. 
| Appendix B.2[ ) as 

H' p « £ Af q , A (a q , A + a^) ]T <"'" r> r (A.12) 
q.A j 
were we have introduced the particle-phonon coupling M qjA given by 



Mq - A = /on 7f E e * q ^p( R ")- (A.13) 

Here /? q denotes the Fourier transform of the modulus square of the Wannier function 

/3 q = J drK(r)|Vi r . 

The function /3 q accounts for the localization of the particles with a hnite width in the 
static trapping potential and thus, for small q, approaches 1. Similarly, for the class of 
potentials we consider (and discuss) in Section^ the Fourier transform of the particle- 
crystal potentials approaches a hnite value at q = 0. Therefore the particle-phonon 
coupling behaves like q 1 / 2 for small momenta, while it is large for q ~ n/a. 

The Hamiltonian of the entire model is given by iJ e ff = H c + H p + H' cp and reads 

H ^ = - J E C J C J + \ E V iA c ] c 3°i 
(id) »»i 

+ E M q^e lqR? c]c J (a q , A + a t _ q , A ) + ^^ q ,Aa qiA a q a. (A.14) 
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The crystal motion given by Hamiltonian (j A. 14|) corresponds to a set of uncoupled 
harmonic oscillators under the influence of an extra particle density dependent force. 
Similarly to the problem of a charged harmonic oscillator in a constant electric field [58| 
such a force displaces the crystal molecules from their original position proportional 
to its strength while keeping the oscillation frequency fixed. A crystal molecule at 
position Ri is therefore displaced by 

7 ^«'" 'A, (A.15) 

with the Fourier transform of extra particle density ?/ q = e~ iqr -»ctcj. 

If the relative displacement between two crystal molecules <5vy = Vj — Vj at 
neighbouring sites i,j is small on the scale of the lattice constant (#v.y) <C a this 
effect can be neglected. 

Appendix B. Phonon spectrum and particle-phonon coupling 

In this section we derive the specific form of the dispersion relation for one- and 
two-dimensional dipolar crystals as well as the crystal-particle interaction. 

Appendix B.l. The phonon spectrum 

The crystalline phase is characterized by small displacements of the molecules from 
their equilibrium positions R^ = R^ + u, . A dipolar crystal that is trapped in one or 
two dimensions by an optical trap with a trapping frequency loj_ is described by the 
following Hamiltonian 

^=ES + Sf^ + ?t5A# (B1) 

where denote the momenta of the molecules. The dispersion relation is found 
from the second order correction of the expansion of the molecule- molecule interaction 
potential 

vW(R i -R j ) = D ij ( Ui -u j ) 2 (B.2) 

with Da = W CC (R°-R°). (B.3) 

z 

We make an ansatz for the displacement as 



»< = ^m^r/^°^ + <'-^ (B ' 4 > 

q,A V M ' 

and can identify a q x and a qi A as creation and annihilation operators of phonons with 
the polarization e\ and quasi momentum q iff the matrix Dij is diagonal. Then the 
polarization vectors for phonons with a dispersion ujq t \ are the eigenvectors of Dij. 
With the discrete Fourier transform 

q q 
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we can write Hamiltonian (jB.ip as 

Hc = 2^rH( p q p -q + m c w q u q u -q) = H^q,A ( a lx a ^ + l) 

C q q,A V 1 

and find the dispersion relation from calculating the eigenvalues of D^. The +1/2 
contribution in the last equation is the vacuum zero point energy, a constant energy 
shift that is ommitted henceforth. 



Appendix B.l.l. Phonon spectrum in ID: In a ID crystal tube the second order of 
the expansion in the displacement of the interaction potential is given by 

r ( ( 2! (R, - R,-) = <f; I () 3 i, I , ( "T^£ . (B.6) 




In g-space we find 



R R 0|5 ^„-,"v" ,• (B- 7 ) 



3 



with / 9 = 5j4sin 2 (gaj/2)/j 5 = 2C(5)-iz5( e ' i9Q )-^ 5 (e^ a ), (B.8) 

j>0 

with C(5) the zeta function at 5 and Li§(x) the polylogarithm of fifth order at x. Here 
we have used the fact that in the crystal tube the equilibrium positions are given by 
R? = ai with the lattice spacing a. 

Since Dij is diagonal in the canonical basis the polarization vectors are given 
by the canonical basis vectors and when including the optical trapping potential, the 
second term on the right handside of equation (|B.1[) , the dispersion relations are found 
as 

/ 1 1 A2 t 

(B.9) 



u q ,± = \ UJ 2 X - (B.10) 

The longitudinal dispersion relation u) q ii is acoustic while the two transversal ones 
denoted by uj q .± show an optical behaviour. The fact that the transversal dispersion 
relations can become imaginary points towards an instability of the crystal, reflected 
by the requirement for a strong transversal trapping, see Section [3l Equations (|B.9|1 
and (jB.lOp show that the optical (transversal) modes decouple from the longitudinal 
one if 




is fulfilled. 



/I5dg 
a 5 m c 
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Appendix B.1.2. Phonon spectrum in 2D: We consider a dipolar crystal that forms 
in the xy-plane. The expansion of the molecule molecule potential to second order in 
the displacement gives 



^cc - | AR0 |7 



5AR9 jiX AR9 jiy 5AR°J - |AR?/ 



V |AR°, 



_ DO 



where we have used ARy = R^ — R , . AR, : _ = R" — R" and similarly ARij X = 
R° — R?j iX , ARij :V — R® y — Rj t y introduce i . = i — j and can write It is difficult to 
diagonalizc the matrix in the last equation. We can however use the block diagonal 
form of the matrix in the last equation and, including the optical trapping potential 
as we have done in the ID case above, write the dispersion relations as 




(B.12) 

y u~//i c 

where denote the two eigenvalues of the xy-block of the matrix and f q is given 
by equation (|B.8|) . The eigenvalues give a longitudinal and a transversal acoustic 
dispersion relation. As noted above we cannot write them in a closed form but an 
approximate result may be obtained by including only particles in the interaction 
that are within some finite range of each other. We show a numeric evaluation of the 
dispersion relation in figure [3J 



Appendix B.2. The particle-phonon coupling 

We denote the crystal-particle interaction potential by V cp (r — Uj) where r and Rj 
denote the position of an extra particle and a crystal particle respectively. The full 
potential is expanded up to first order in the displacement R., = R° + Uj as 

V cp (r)= £ Uj -W cp (r-R°) (B.13) 

3 

= T^Ad E fike i(k - q)r V q %(q) (B.14) 

where d denotes the dimension of the setup and we have used the discrete Fourier 
transform of the interaction potential V r cp (R°) = ^ q e lqR; 14p(q)/v / 27r and the 

displacement uj — ^ k e jkR i u\ L j\f 7 2jx . The Hamiltonian is found by integration of 
the extra particle density p(r) over the interaction H\ = J drp(r)V cp (r) which gives 

H 1 = V2^ d ^p(q)T/ cp (q)qu q . (B.15) 
q 

The extra particle density p(r) is defined through the single particle operator, which 
reads in site representation ip( r ) — Em c m ( f>m( r ) with </> m (r) the wavefunction of an 
extra particle at site m, and the annihilation operator c m as p(r) = ip* (r)ip(r). In the 
tight binding limit the particles are strongly localized at sites and the overlap of the 
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wavefunctions of particles at different sites is neglected. This approximation gives 

p(q) « ^4c m / -^L|0 m (r)|V<* r (B.16) 

= E^c™ / ^|0o(r)|V^> (B.17) 



'2-K 



^E^^- (B.18) 



We can work out /3„ = f G?r|0(r)| 2 e ,qr by making a separation ansatz for the 
wavefunction </>(r) = (fi x (x)(j)y(y)(j) z (z) and thus /3 q = /3q x Pq v /3 qz . In the directions 
of confinement the wavefunction is taken to be Gaussian 4> z (z) = e~ z ' ! / 2a i / t -K 1 ^o}[ 1 
with a± = ^h/m p u±_ which gives (i qz — e~ a i' 3 ^ 4 . 4>a{r) is taken a Wannier function 
in all other directions thus /3 q = J dr|uJo(r)| 2 e lqr with wo(r) the Wannier function of 
the lowest Bloch band. The interaction Hamiltonian is 

#1= EE 4c m e^/3 q V;p(q)qu q (B.19) 

q m 

= E ^q,Ae' qr ™4 Cm (a q , A + at q>A ) (B.20) 

m, q ,A 

where we can use (|A.7|) to identify 

M q , A - (? l /2iVm c c I ; q ^) 1 /2 / 3 ql / cp ( q ) q e A . (B.21) 

Extra particles interacting with crystal molecules will force the latter to new 
equilibrium positions displaced from their original positions by an extra particle 
density dependent displacement 



v = 2 y 1 e <q(H5-rS) c t (B . 22) 

The impact of the displacement is neglected if the relative shift between two 
neighbouring molecules |v$ — Vj+^l = (vj — vt+s) is small compared to the lattice 
constant. 

Appendix C. Correlation functions, corrections and diagonalization of the 
master equation 

Appendix C.l. Expectation values and correlation functions 

The crystal is in a thermal equilibrium at temperature T with a reference state given 

by 



Pb 

q,A 



b = jjexp [-^r« q .A a q,Aj (1 - exp 



hie, 



q,A 



With the displacement of crystal molecules by an extra particle located at site r ! 



3 ' 

X, = exp[-E^,Ae iqr V_ q)A ~ «q,A)], (C.l) 

q,A 
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we find 

x\Xt = exp [5> q ,A(e iqr ° fc ~ ^ r ?)(a q , A - al q>A )] 

q,A 

= exp[^(Z^a q ,A-^ A *«q,A)] =\{D{a^Z^\ (C.2) 

q.A q,A 

where we have introduced Z^ x = u q . A (e~ 4qr < 1 — e~ 4qr °) and the displacement operator 
D(a,a) = e Qat ~" a . In Section \2.'2\ we have introduced the (thermal) expectation 
value of a bath operator O by (O) = ^B,{Op B } where tre denotes the trace over the 
bath degrees of freedom. The expectation value of the displacement operator is given 
by (£>(a q , A ,a)) = exp[-(n q , A (T) + ±)|a| 2 ] where n q , A (T) = l/(exp[/iw q , A /feT] - 1) 
denotes the thermal expectation value of the phonon number operator. Therefore we 
find 

{X£x t ) = e- ST (C.3) 

with 



S T = 2^< A sin [|(r° - r°)J (2n q , A (T) + l) 

q,A 



(C.4) 



Notice that the link — is always a symmetry axis of the integration interval, the 
first Brillouin zone. Since itq iA and n qjA (T) are invariant under a rotation with the 
symmetry of the Brillouin zone the sum over q in equation (|C.4[) becomes independent 
of the orientation of that link. 

Expectation values of time dependent displacements can we worked out in 
a similar fashion. With the explicit form of the bath Hamiltonian H-q = 
Ylq x ^q.A&q A a q.A we can write a bath operator in the interaction picture as a,q y \(t) = 
<2q, A e IWc| A * and wc find 

xl{t) Xl {t) = exp [J2 (^' A « q .A« - Z^i X *a^{t)) 



q,A 



= \{D{a^,Z^e-^). 

q,A 



(C.5) 



In the following we introduce Z?f(t) = Zffe 1 ^ 1 and use the identity for 
displacement operators P D(aq, a)D(a p , (3) = Yl a/2)Z?(a q , a + j3) to 

calculate 

(X^X^X^t - T)Mt - r)) =Ul[(D(a^ x , Z%\t))D(a p , K , 2%>(t - r))) 

q,A p 7 K, 

= \{{D{Z^ X {t - r), (i)/2)D(a q , A , Z%\t) + Z%\t - r))) 



q,A 



-2S, 



ex p { - E [K.>( T ) + i)z% x *W% x (t r) 

+ n^(T)Z% X (t)Z^*(t-T)]} 



q.A 



= e~ 2ST e _ *y (x ' T) . 
Here $^-(r, T) is found by inserting for Z(t) as 



(C.6) 



$*?(r, T) = £ <a [K,A + l)9^e-^ + % A <$V<^] (C.7) 



q,A 
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with g$ = (e" lqr ? -e~ l ^){e^ r ° -e lqr *). The quantity $g(r, T) depends only on the 
relative time r = t — t' and the two links rf — r° and r° k — r° . A symmetry argument 
that relies on the fact that (qr°) is a projection on a symmetry axis of the integration 
interval, the first Brillouin zone, while everything else is an even function in q allows 
us to write (IC.7I) as 



(r, T) = <x9% [ coth (^) cobKat) - i sin( Wq , A r)] (C.8) 

'I A 



with ff« = cos[q(r° - r°)] - cos[q(r° - r°)] 

-co S [q(r°-r°)]+cos[q(r°-r^]. (C.9) 

From this form one can immediately read off the important relations <&^- (t, T) = 

$y ; (r,T), $g(-T,T) = ^(rj 1 ) and |*g(T,T)| < 2S T - Equation JUU) is in 
accordance with the literature [58| . 

For many sites we take the continuum limit and replace the summation over q by 
an integration over the first Brillouin zone of volume Vbz- A variable transformation 
where q x is replaced by f(q d ~ 1 ,w), a function of the remaining momenta <7 d_1 and 
w = w q ,A, allows us to introduce the spectral density 



A ^ 



g^(i d ~ 1 : W ) 



$$(t,T) = / d«;J#(«>) coth ( T71 - 7F ) cos(wt) - isinfar) . (CIO) 



where the integration interval depends strongly on w and the form of the Brillouin 
zone while w ranges over all frequencies. In terms of the spectral density equation (IC.8|) 
reads 

/ hw 

The bath correlation functions £y (r, T) that enter the master equation are thus given 

by 

$(r,T) = (xj(t)X j (t)xl(l/)X l (t'))-e- 2S - (C.ll) 

= e -2S T(e -^](r,T)_ 1) (C12) 



and satisfy Q(-r,T) = ^*(r,T) 



Appendix C.2. Corrections to the Master Equation 

The corrections to the master equation are given by (c.f. Section \2~2 



r ii( T )=^l drRe[$(r,T)], (C.13) 



AS(T) 



fi 2 

r2 



/>OC 

/ drRe[$(r,T)], 

JO 

/>OC 

/ drIm[^(r,T)], (C.14) 
Jo 

^( T )=tfJ o dTTReKg(r,T)]. (C.15) 

t2 /-oo 

^( T )=Jfj o dTTlm[$(T,T)], (C.16) 

We can find explicit approximate results to these corrections in the two limiting cases 
of a weak, St <C 1, and a strong, St >■ 1, coupling. 
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(b) 

1.0 




- Regj£(T,T)] 









10 

TlO D 



15 
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Figure CI. Real part of (a) the function 3?S(t, T) and of (b) the correlation 
function £^?(r, T) for a ID model with Af q c< q/ y/ /UJ^ and u; q = l>d sin(g/2) as 
a function of the (dimcnsionless) time tu>y>- Shown are the respective functions 
for the swapping of a particle, i = I and k = j (solid lines), and the hopping 
of particles into the same direction, j = k (dashed lines). Notice the overall 
minimum of the functions \&M[t,T)\ is attained at r = for the swap process, 
i.e. i = I and k = j. In panel (b) we chose St =3/4, while for larger values of St 
the tail of Re[^J(r, T)], together with all other processes, is strongly suppressed. 
In the strong coupling limit (St 3> 1) the latter become of the order of 0(e~ 2ST ). 



Appendix C.2.1. Strong coupling limit The real part of $fj(r, T) has a global 

minimum at r = To, see figure [CTT a). Therefore the function e - * 4 ^ 1 "'^ is strongly 
peaked at this minimum, see figure [ClT bV For a sufficiently strong coupling, St ^ 1, 
all contributions to the r integration come from a close range around the minimum, 
r — to. Then we may approximate the integral by performing a stationary phase 
approximation which features an expansion of (t, T) around tq . We find 

Re{*g(r, T)} = ]T J$(w) coth ( J^L) ( cos( w (t - t )) cos( W t) 

- sin(w(T - t )) sin(Wo)) , 

~ J ij( w ) coth (^) ( cos(wto) - w(t - t ) sin(wT ) 



^w 2 (t - t ) 2 cos(wto)), 



and in the same way 

Im{$g(T, T)} » 5^ coth (gjr^;) ( sin(wTo) - ib(t - t ) cos(u;to) 



iw 2 (T - t ) 2 sm(u;To)). 
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This gives for real and imaginary part of the correlation function 

J.f/M coth (^y) ( cos(u;ro) 

hw 



exp 



[: 



1 



x exp 



[J2 J « M coth (jjjfc^f ) W ( T ~ sin ( WT °) + l w2 ( T - T °) 2 cos(wr ) 



Tn [ 51 J ^ ^ ^ sin(wro) + w(t - t ) cos(u;t )) 



sin 



where real and imaginary part differ only by the trigonometric function at the 
beginning of the last line in the last equation. The first line on the right hand side 
of the last equation does not depend on r. It exponentially suppresses the rest of the 
function by 

hw 



J2 J %( w ) coth 



2k B T 



i(wT Q ) - 1). 



This factor is of the order of St unless tq < 1. In the strong coupling limit where 
St S> 1 every process characterized by J^j(w) becomes therefore strongly suppressed 
if the corresponding correlation is not peaked around r = 0. An analysis of the 
correlation functions shows that this condition is only met by correlations to the swap 
process, i = I and k = j. For these processes real and imaginary part of the correlation 
functions are given by 



Re 
Im 



[$(r,T)] « exp (^J#H coth 



hw 

2k B T 



2 2 

W T 



Tn [E-#("H< ai7 > 



and the corrections (|C.13[) - (|C.16|) are simply Gaussian integrals that calculate as 



F io 
1 01 



(T) 



7o 1 ?(T) 



A™(T) 



J2 n l/2 e -B*/4A T 

fi 2 ~2 y/Ar~ 

j2 n l/2 e ~B 2 /4A T 



h 2 



Erfi(S/2 v / Ir), 




1 Be- s2 / 4AT V^Erfi( J B/2VA^) 



2A T 

J 2 7rV2 Se- s2 / 4A ^ 



4A 



3/2 



(C.18) 
(C.19) 
(C.20) 



ft 2 4 



.4 



3/2 



(C.21) 



where Erfi denotes the Error function and the quantities At and B are given by 



B 



1 



dw-J^(w)w coth 



g?u> Jq°(u>)u>. 



/ hw 
\2k B T 



(C.22) 
(C.23) 



In the strong coupling limit we can do a simple estimate to show how AJJ(T) and 
Y01CO are part of an J/E p expansion for ID systems. The expansion parameter 
becomes visible when including only first order in the series expansion of — $>q1(t,T) 
around its maximum. This is equivalent to the limit of At going to zero. In this limit 
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we find 

->°° m - (C - 25) 

The quantity B for this process reads 

B = J dwJ™(w)w = J dq(^ff sin(g/2) 3 (C.26) 

~ -T— / dff^tl " cos(g)] = M n t M+1 . (C.27) 

In the two models in Section 2] we find that the onsite phonon mediated interaction 
— 2E p dominates over all the others. Therefore we can approximately write 
B w Ep/nhu>D with which we find 

Aj°(T)«7r^, (C.28) 
E p 

7o?(r) «^ 2 ^- (C29) 

Appendix C.2.2. Weak coupling limit In the weak coupling limit, St <C 1, we use the 
identity (t, T) < 2St to perform a pointwise expansion of the correlation functions 
(ICT21) in $§(t,T) as 

4?(r,T)«e- 2S -^J(r,T). (C.30) 

With this expansion real and imaginary part of (r, T) simply become 

Re[^(r,T)] = e" 25 - | d«;J#(«;) coth {^^j cos(wr), (CM) 
Im{$ (t,T)} = - er 2ST f dwJ^Hw) sm(wr), (C.32) 



where the w integration ranges over all frequencies. This way the Debye frequency of 
our model functions as a natural cutoff (compare [68]). The functional identities 



/ drcos(wr) = irS(w), 
Jo 

POO 

I drsin(wr) = P(l/w), 
Jo 

/>oo 

/ drr cos(wr) = —nS'(w), 
Jo 

/>00 

/ drTsm(wT) = - P(l/w 2 ), 
Jo 



with P(x) the Cauchy principal value of x and the delta functional 5(x), allow us to 
find the corrections in the weak coupling limit as 

^(M)-^[^Mcoth(^)]_ o , (C.33) 
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Ag(r,t) - ^ lim| dwJ%(w)-^, (C.34) 
7#(r,t) = - ^lim/ ^Ncoth(^) I |^ I (C.35) 



$(r,t) = - tjtt^^HI^q, (C.36) 



h 2 ' 

where we have written out the Cauchy principal value integrals. We remark that the 
integration (|C.35|) becomes infinite in the limit of zero temperature. 



Appendix C.3. Diagonalization of the master equation 

As noted above the corrections to the master equation Iy (T), (T), jfJ(T) and 
5fj(T) are just matrix entries and have little physical meaning by themselves. 

In this section we show how to estimate the energy of the corrections to the master 
equation (fi"6)) for a single particle at many sites. We start by treating the self energies, 

£ AgCTJMw = EE A:^ +m+ "4c l+m+ „ (C.37) 

(ij){kl) i m,n 

= J2J2 A™ m ™+V<^ +r ° >ct Cq (C.38) 

q m,n 

where the sums over m, n range over basis vectors in the lattice. In the second step we 
have made use of the fact that the energies A t i + i " l ^ l +m+n do not depend on the index 

i and that in q space we can write X)i c l c «+™ = E q e ' qr "' c q c< )' ^ e e ig enva l ues to 
(|C.38|) are simply given by 

A q = A^ + V q(r ™ +r ° ) (C.39) 

and the largest eigenvalue gives an upper bound to the energy of this term in the 
master equation. 

To determine the energy of the thermal fluctuations we notice that this term can 
be written as 

i ]T ta?V(T)({b ij b kl ,p s (t)}-2b kiPs (t)b ij ) (C.40) 

(ij)(kl) 

where we have introduced bij = c\cj and we estimate its amplitude by the energy of 
J2{ij)(ki) f^ij(T)bijbki- Since T^j (T) has the same properties as A*j- (T) the eigenvalues 
are given by 

r q = E r™r + "e 4q(r '" +r " ) - (C41) 

in, 7i 

The dissipative term proportional to j^j(T) in the single particle limit is given by 



E 7^(T)[6 4 ,(E^-E^') PsW+Ps(<) (^ &fc ''~^^')^ 

v V fc' v 

-^■psw(E^'' -E 5fc - (E 5fe '' -E^'W*)M' ( c - 42 ) 



h 

iij)(hi) 
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where k' and I' denote the nearest neighbours of k and I. We can write the first two 
terms as 



[ij)(kl) k' V i m,n,o 

i+m,i-\-m-\-n-\-o /rri\ \ t / r\ A o\ 

-id+m (T)Jclc l+m+n+0 , (C.43) 

E 7?j (?) ( E ^ ^ E bkl ') b V = E E (lk+m+n^+m+n+o 



(ij)(kl) k' I' k m,n,o 



k,k+m \ t (r< aa\ 

lk+m+n,k+m+n+o J c k c k+rn+n+o, {^.44} 



and as before the sums over m, n, o range over all basis vectors. The property 
lij (T) — J^i (T) allows us to collect the first two terms in (IC.42[) by an anticommutator 

~ Y\ 2^ \Y0,m ( T )-7o,m (T)J }^ C l C t+rr l +n+o,Ps(t)j- (C.45) 

m,n,o i 

In this case the amplitude of the entire term can be estimated by the eigenvalues of 
the operator inside of the anti commutator. They are given by 

7 q = E {loX n ' m+n+0 ( T )-^ n+ °^ (C46) 

m,n,o 

We use the same line of argumentation to estimate the eigenvalues of the dissipative 
correction proportional to 5fj and find 

<*q= E (<^+" ,ro+n+ °(T) - 5^ (C47) 

m,n,o 

In a one-dimensional system the extra particles locate at sites of postion r® m = ma 
with the lattice constant a. There are only the two basisvectors ±a and therefore the 
indices m,n and o range only over ±1. This gives 

r,(T) = r™^ + "(r)e jq ( r ™ +r ") (c.48) 
= rJ;?(T)e a «° + rJ;?(T) + r^^(T) + r a ^ 2 (T) e ^ ,;2 « a (c.49) 

= 2Tj;?(T) + 2rJ; 2 (T) cos(ga). (C.50) 

In the second step we have used the fact that the origin of the indices i,j, k, I is the 
function gfj(q) and therefore the corrections show the same symmetries with respect 
to the indices as g^iq) itself. In the same way we find the other eigenvalues as 

A,(T) = 2[Aj?(T) + AH(T) cos(qa)}, (C.51) 
7 ,(T) = 2[( 7o 2 j 5 (T) - 7 $(T)) cos(3 g a) 

+ (lolm + 7o ° 1 1 (T) + 7o V 1 (T) - ToV'CO) cosM], (C.52) 
S q (T) =2[(6l 3 1 (T)~S 1 3 1 (T))co S (3qa) 

+ (S 2 ol(T) + S° Q l(T) + d° Q f\T) - Slf^T)) cosM]. (C.53) 
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